﻿/*说明：
本方法仅作为牛顿迭代之外的另一种尝试，之所以又写这个算法，是因为发现numpy、MATLAB等主流的求解工具在求解多项式方程时
入手的方向都是将方程的伴随矩阵求出来，进而利用其他方法将伴随矩阵的特征值求出，而求出的特征值即为多项式方程的根。
本方法是在初赛的结果已经完全满足比赛要求之后才开始的一个新探索，仅供参考。
结果：
1.本方法在求解重根时的效果依然不佳，误差常常高达到1e-2。
2.本方法在引入高精度库之后仍然没法将误差下降至指定要求。
3.本方法相较于牛顿迭代的优势在于其能够求出方程的复平面的根。

注：numpy和MATLAB在求解多项式遇到重根时的表现同样是不佳，与本方法的误差相差不大。
*/

// 自定义头文件
#include "template_problem01_characristic.hxx"

// ACIS
#include "acis/acis.hxx"
#include "acis/api.err"
#include "acis/api.hxx"
#include "acis/body.hxx"
#include "acis/check.hxx"
#include "acis/cstrapi.hxx"
#include "acis/logical.h"
#include "acis/module.hxx"
#include "acis/primtive.hxx"

#include "acis/law_base.hxx"
#include "acis/law.hxx"
#include "acis_utils.hpp"
#include <iostream>
#include <complex>
#include <vector>
#include <cmath>
#include <algorithm>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector_expression.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/math/ccmath/floor.hpp>
#include <boost/numeric/ublas/io.hpp>//用来输出矩阵

using namespace boost::numeric::ublas;
using namespace boost::multiprecision;
using namespace boost::math::ccmath;
//注意，在本文件中的所有vector都是在boost::multiprecision命名空间中的vector，
//而不是std中的vector

const cpp_dec_float_50 epsilon = 1e-10;
const int MAX_ITER = 1000;


/**
 * @brief		获取方程的伴随矩阵
 * @param[in]	coef	方程系数数组
 * @param[in]	N		方程最大次数
 * @return		返回方程的伴随矩阵
 */
matrix<cpp_dec_float_50> buildCompanionMatrix(const std::vector<cpp_dec_float_50>& coef, int N);

/**
 * @brief		用Gram-Schimidt正交法进行QR分解
 * @param[in]	A	方程的伴随矩阵
 * @param[out]	Q		Q矩阵（正交）
 * @param[out]	R		R矩阵（上三角）
 */
void qrDecomposition_GS(const matrix<cpp_dec_float_50>& A, matrix<cpp_dec_float_50>& Q, matrix<cpp_dec_float_50>& R, int n);

/**
 * @brief		用Householder转换进行QR分解
 * @param[in]	A	方程的伴随矩阵
 * @param[out]	Q		Q矩阵（正交）
 * @param[out]	R		R矩阵（上三角）
 */
void qrDecomposition_H(const matrix<cpp_dec_float_50>& A, matrix<cpp_dec_float_50>& Q, matrix<cpp_dec_float_50>& R, int n);

/**
 * @brief		用QR方法将根迭代出来
 * @param[in]	A	       方程的伴随矩阵
 * @param[out]  eigenvalues求出的特征值，即为方程的根
 * @param[in]	max_iter   最大迭代次数
 * @return		方程的所有根
 */
void qrAlgorithm(matrix<cpp_dec_float_50>& A, std::vector<std::complex<cpp_dec_float_50>>& eigenvalues, int n);




// QR分解——Gram-Schmidt正交方法
void qrDecomposition_GS(const matrix<cpp_dec_float_50>& A, matrix<cpp_dec_float_50>& Q, matrix<cpp_dec_float_50>& R, int n) {
	Q = matrix<cpp_dec_float_50>(n, n);
	R = matrix<cpp_dec_float_50>(n, n);


	//习惯用i代表行，j代表列
	for (int j = 0; j < n; j++) {
		std::vector<cpp_dec_float_50> v(n, 0.0);
		for (int i = 0; i < n; i++) {
			v[i] = A(i, j);
		}

		for (int k = 0; k < j; k++) {
			cpp_dec_float_50 dot = 0.0;
			for (int i = 0; i < n; i++) {
				dot += Q(i, k) * v[i];
			}
			for (int i = 0; i < n; i++) {
				v[i] -= dot * Q(i, k);
			}
		}

		cpp_dec_float_50 norm = 0.0;
		for (int i = 0; i < n; i++) {
			norm += v[i] * v[i];
		}
		norm = sqrt(norm);

		for (int i = 0; i < n; i++) {
			Q(i, j) = v[i] / norm;
		}


		for (int i = 0; i <= j; i++) {
			R(i, j) = 0.0;
			for (int k = 0; k < n; k++) {
				R(i, j) += Q(k, i) * A(k, j);
			}
		}

	}

}

// QR分解——Householder正交方法
void qrDecomposition_H(const matrix<cpp_dec_float_50>& A, matrix<cpp_dec_float_50>& Q, matrix<cpp_dec_float_50>& R, int n) {
	//Q为单位阵
	Q = identity_matrix<cpp_dec_float_50>(n);
	R = A;

	for (int i = 0; i < n - 1; i++) {

		vector<cpp_dec_float_50> x(n - i);
		for (int j = i; j < n; j++) {
			x[j - i] = R(j, i);
		}
		vector<cpp_dec_float_50> e(n - i);
		if (x[0] >= 0) {
			e[0] = norm_2(x);
		}
		else {
			e[0] = -norm_2(x);
		}
		vector<cpp_dec_float_50> u = x - e;
		vector<cpp_dec_float_50> v(n - i);
		cpp_dec_float_50 u_norm = norm_2(u);
		for (int j = 0; j < u.size(); j++) {
			v[j] = u[j] / u_norm;
		}

		matrix<cpp_dec_float_50> H = identity_matrix<cpp_dec_float_50>(n);
		project(H, range(i, n), range(i, n)) = project(H, range(i, n), range(i, n)) - 2 * outer_prod(v, v);

		// 更新 R
		R = prod(H, R);

		// 更新 Q
		Q = prod(Q, trans(H));

	}

}


// QR方法求特征值
void qrAlgorithm(matrix<cpp_dec_float_50>& A, std::vector<std::complex<cpp_dec_float_50>>& eigenvalues, int n) {
	//int n = A.size();
	matrix<cpp_dec_float_50> Q, R;
	matrix<cpp_dec_float_50> E(n, n);
	cpp_dec_float_50 delta = 0;
	vector<std::complex<cpp_dec_float_50>> roots(n);//用于保存最终的根
	vector<std::complex<cpp_dec_float_50>> temp_roots(n);//用于收敛判断

	for (int iter = 0; iter < MAX_ITER; iter++) {

		delta = A(n - 1, n - 1);
		for (int i = 0; i < n; i++) {
			E(i, i) = delta;
		}
		A = A - E;

		qrDecomposition_H(A, Q, R, n);
		//matrix<cpp_dec_float_50> tempA = A;

		// 更新A = R * Q
		//A = prod(R, Q);
		A = prod(R, Q) + E; //shift方法

		// 提取特征值
		temp_roots = roots;
		roots = vector<std::complex<cpp_dec_float_50>>(n);
		int i = 0;
		while (i < n) {
			if (i<n - 1 && abs(A(i + 1, i))>epsilon) {
				//检测到2*2子块,准备解二次方程
				cpp_dec_float_50 a = 1;
				cpp_dec_float_50 b = -A(i, i) - A(i + 1, i + 1);
				cpp_dec_float_50 c = A(i, i) * A(i + 1, i + 1) - A(i + 1, i) * A(i, i + 1);
				std::complex<cpp_dec_float_50> root1;
				std::complex<cpp_dec_float_50> root2;
				std::cout << "判别式为：" << b * b - 4 * a * c << std::endl;
				if (b * b - 4 * a * c > 0) {
					root1 = std::complex<cpp_dec_float_50>((-b + sqrt(b * b - 4 * a * c)) / (2 * a), 0);
					root2 = std::complex<cpp_dec_float_50>((-b - sqrt(b * b - 4 * a * c)) / (2 * a), 0);
				}
				else if (b * b - 4 * a * c < 0) {

					root1 = std::complex<cpp_dec_float_50>(-b / (2 * a), sqrt(4 * a * c - b * b) / (2 * a));

					std::cout << sqrt(4 * a * c - b * b) / (2 * a) << std::endl;
					std::cout << root1 << std::endl;
					root2 = std::complex<cpp_dec_float_50>(-b / (2 * a), -sqrt(4 * a * c - b * b) / (2 * a));
				}
				//cpp_dec_float_50 root1 = (-b + sqrt(b * b - 4 * a * c)) / 2 * a;
				//cpp_dec_float_50 root2 = (-b - sqrt(b * b - 4 * a * c)) / 2 * a;
				std::cout << "系数a,b,c分别为：" << std::endl << a << std::endl << b << std::endl << c << std::endl;
				std::cout << "结果分别为：" << std::endl << root1 << std::endl << root2 << std::endl << std::endl;
				roots[i] = root1;
				i = i + 1;
				roots[i] = root2;
				i = i + 1;
			}
			else {
				roots[i] = std::complex<cpp_dec_float_50>(A(i, i), 0);
				i = i + 1;
			}
		}

		// 检查是否收敛
		// 用根的变化来判别是否收敛
		//cpp_dec_float_50 sum = 0;
		//for (int j = 0; j < n; j++) {
		//	sum += abs(temp_roots[j] - roots[j]);
		//}
		//if (sum < epsilon) {
		//	std::cout << iter << std::endl;
		//	break;
		//}

		//shift方法的收敛条件
		cpp_dec_float_50 sum = 0;
		for (int i = 0; i < n - 1; i++) {
			sum += A(i + 1, i);
		}
		if (abs(sum) < epsilon) {
			std::cout << iter << std::endl;
			break;
		}

		////Frobenius 范数
		//matrix<cpp_dec_float_50> mineA = tempA - A;
		//cpp_dec_float_50 s = 0;
		//for (int i = 0; i < n; i++) {
		//	for (int j = 0; j < n; j++) {
		//		s += mineA(i, j) * mineA(i, j);
		//	}
		//}
		////std::cout << s << std::endl;
		//if (s < epsilon) {
		//	std::cout << iter << std::endl;
		//	break;
		//}

	}

	std::cout << "QR方法得到的矩阵为：" << std::endl;
	for (int j = 0; j < n; j++) {
		std::cout << row(A, j) << std::endl;
	}
	std::cout << "最终求出的根为：" << roots << std::endl;

	for (int i = 0; i < n; i++) {
		eigenvalues.push_back(roots[i]);
	}

}

//获取伴随矩阵
matrix<cpp_dec_float_50> buildCompanionMatrix(const std::vector<cpp_dec_float_50>& coefficients, int n) {
	//Boost.uBLAS库在构造矩阵时默认为零矩阵
	matrix<cpp_dec_float_50> companion(n, n);

	project(companion, range(1, n), range(0, n - 1)) = identity_matrix<cpp_dec_float_50>(n - 1);

	for (int i = 0; i < n; ++i) {
		companion(0, i) = -coefficients[n - 1 - i];
	}

	return companion;
}

outcome answer_api_nroots_of_chara(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int* size, cpp_dec_float_50** answer, int N) {
	API_BEGIN


		//由于矩阵的算法结果为复数，因此不参与测试，这里提供方程的系数用以验算，从上往下次数依次增加
		std::vector<cpp_dec_float_50> vec = { 938575.71754547941070516676048784491090759366872078,
											-149487.85793239880204977344815597956273535645413315,
											7276.4867389446535339778607195391173397941869876344,
											-143.17811216841479104355104462713921918200318488595,
											1 };


	// 构造伴随矩阵
	matrix<cpp_dec_float_50> companion = buildCompanionMatrix(vec, N);

	std::cout << "Companion Matrix: " << std::endl << companion << std::endl;
	//printMatrix(companion);

	// 使用QR方法求特征值
	std::vector<std::complex<cpp_dec_float_50>> roots;
	qrAlgorithm(companion, roots, N);

	std::cout << "roots of the polynomial: " << std::endl;
	for (const std::complex< cpp_dec_float_50>& root : roots) {
		std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << root << std::endl;
	}
	std::cout << std::endl;




	API_END
		return result;
}